Glassy Dynamics in the Adaptive Immune Response Prevents Autoimmune Disease 
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The immune system normally protects the human host against death by infection. However, 
when an immune response is mistakenly directed at self antigens, autoimmune disease can occur. 
We describe a model of protein evolution to simulate the dynamics of the adaptive immune response 
to antigens. Computer simulations of the dynamics of antibody evolution show that different evolu- 
tionary mechanisms, namely gene segment swapping and point mutation, lead to different evolved 
antibody binding affinities. Although a combination of gene segment swapping and point mutation 
can yield a greater affinity to a specific antigen than point mutation alone, the antibodies so evolved 
are highly cross-reactive and would cause autoimmune disease, and this is not the chosen dynam- 
ics of the immune system. We suggest that in the immune system a balance has evolved between 
binding affinity and specificity in the mechanism for searching the amino acid sequence space of 
antibodies. 

PACS numbers: 87.10.+e, 87.15.Aa, 87.17.-d, 87.23.Kg 



One of the main functions of the immune system 
is to generate B cells that secrete antibodies, protein 
molecules that recognize and bind to antigen. The im- 
mune system has evolved a hierarchical strategy to pro- 
duce antibodies that combat invading pathogens Q] . The 
first level is VDJ recombination, in which gene fragments 
are rearranged to encode functional antibodies, produc- 
ing a range of antibodies with a diversity on the order 
of 10 12 to 10 14 . Upon invasion of the body by a foreign 
antigen, the second level of B cell somatic hypermutation 
evolution is initiated, in which division, mutation, and 
selection of B cells occurs. Those B cells that produce 
antibodies that bind the antigen with higher affinities are 
selected and propagated. Thus, the individual point mu- 
tations of somatic hypermutation essentially perform an 
optimizing local search of amino acid sequence space. 

When an immune response is mistakenly directed at 
self antigens, autoimmune disease occurs In most 
cases, the antibody/antigen interaction is highly specific. 
Sometimes, however, the antibody evolved in response 
to one antigen can bind other antigens, and this phe- 
nomenon is called cross- reactivity |2j. Cross- reactivity 
happens when the other antigen has chemical features in 
common with the original antigen and is quantified ex- 
perimentally by measuring the affinity of the antibody 
for the other antigen. Cross-reactivity is one mechanism 
by which autoimmune disease may develop. 

In this Letter, we present two major results. First, 
we show how different evolutionary mechanisms influ- 
ence the relaxation dynamics of an evolving population 
of proteins. Specifically we study the dynamics of an- 
tibody generation by the immune system in response to 
pathogen invasion. These dynamics are crucial to the effi- 
cacy of the immune response. The hierarchical structure 
of our model protein Hamiltonian plays a critical role 
in the dynamics. The hierarchical structure also distin- 
guishes our model from traditional short- or long-range 



spin glass energy models 0, 0|. Second, we show that 
antibodies with significantly higher affinities to antigen 
than those produced in a primary immune response can 
be found with a biologically-plausible evolution process, 
but are more cross-reactive and would greatly increase 
susceptibility to autoimmune disorders. The concept of 
cross-reactivity is related to that of the chaos exponent 
in traditional spin glasses 0, Q, except that for cross- 
reactivity we are interested in the immediate energy re- 
sponse rather than in the equilibrated response of the 
system to a change in the couplings. Taking our two 
results together, we show that in the primary adaptive 
immune response a careful balance has evolved between 
specificity for and binding affinity to foreign antigen. 

To express the interaction energy between antibodies 
and antigen, we use the generalized NK model, a model 
that captures important features of the immune response 
to vaccination and disease and of protein molecular evo- 
lution 0, S Q- The energy of interaction for a given 
protein or antibody sequence is defined in the model as 
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The term U sd includes interactions within a secondary 
structure. Interactions between secondary structures are 
essential to protein folding and function. These terms are 
in t/ s . The total number of interactions with a typi- 
cal amino acid is roughly 12, and half of these are in U sd 
(2(K-1)), and the other half in U sd ~ sd (D{M-l)/N). 
The term U c is the direct binding interaction between 
the antibody and antigen. The number of antibody sec- 
ondary structural subdomains is M = 10, and the num- 
ber of antibody amino acids involved directly in the bind- 
ing is known to be approximately P = 5 [8j. We define 
the binding constant between antibody and antigen to 
be equal to K = exp a ~ fc ^ where a and b are constants 
found by comparison to immunological data @. Thus, 
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the lower the overall energy, the better the binding affin- 
ity of the antibody to the given antigen. The secondary 
structure energy U sd is 
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u sd = 



y/M{N - K) 



(2) 

where N = 10 is the number of amino acids in a sub- 
domain, and K = 4 [10| specifies the range of the local 
interactions within the secondary structure. We consider 
five chemically distinct amino acid classes (negative, pos- 
itive, polar, hydrophobic, and other), and all subdomains 
belong to one of L = 5 different types (helices, strands, 
loops, turns, and others). The quenched Gaussian ran- 
dom number a ai is different for each value of its argument 
for a given subdomain type, o^. All of the Gaussian a 
values have zero mean and unit variance. The energy of 
interaction between secondary structures i and j is 
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where D = 6 specifies the number of interactions between 
secondary structures. Here cry and the interacting amino 
acids, h, ...,Iki are selected at random for each interac- 
tion k. The chemical binding energy of each antibody 
amino acid to the antigen is given by U? = Oi(ai) / 'y/P '. 
The contributing amino acid, cii, and the unit-normal 
weight of the binding, c^, are chosen at random. 

In our first study, we replicate the immune system 
dynamics. Each protein sequence is of the length of 
N x M — 100 amino acids, corresponding to the length 
of the heavy chain variable region of an antibody. Ini- 
tially there are 10 3 sequences, as in the human immune 
response Two different strategics arc used in our sim- 
ulations to search protein sequence space for high affinity 
antibodies. The first strategy mimics the normal adap- 
tive humoral immune response. It starts with the combi- 
natorial joining of optimized subdomains, after which the 
sequences undergo rounds of point mutation (PM) and 
selection. This is a simulation of the VDJ recombination, 
somatic hypermutation, and clonal selection that occurs 
in B cell development Q- We generate five optimized 
subdomain pools, each composed of 300 low-energy sub- 
domains, corresponding to the L — 5 types. We use three 
sequences from each subdomain pool in our VDJ recom- 
bination, mimicking the known diversity [jj. 

In our second study, we include a more powerful search 
of sequence space in the antibody evolution dynamics. 
VDJ recombination is used to generate an initial popu- 
lation, as in the normal immune response, but in addi- 
tion to somatic hypermutation, we perform gene segment 
swapping (GSS) during each round of evolution and selec- 
tion. GSS-type processes are used by experimentalists to 
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FIG. 1: Evolution of the affinity energy for the cases of point 
mutation (PM) only and gene segment swapping (GSS) plus 
PM as a function of the number of rounds of mutation and 
selection used to evolve the population of antibodies. Energies 
evolved at larger number of rounds shown in inset. 



produce antibodies with binding constants s» 10 11 — 10 13 
1/mol [HI, 

and exist within the natural hierarchy of evo- 
lutionary events 0, ^3 • 

In the process of GSS, a subdomain belonging to one 
of the five types is replaced by another one from the 
optimized subdomain pool of the same type. Each se- 
quence undergoes an average of 0.5 point mutations in 
each round of selection in both strategies. In GSS+PM, 
each subdomain in a sequence has a probability of 0.05 
of being replaced in the process of GSS. Following mu- 
tation, selection occurs, and the 20% lowest energy anti- 
bodies are kept and are amplified to form the population 
of 10 3 sequences for the next round of mutation and se- 
lection. The primary response is comprised of 30 rounds 
of affinity maturation, corresponding to a lag phase of 
approximately 10 days [t| , during which B cells undergo 
clonal selection in response to antigen and differentiate 
into plasma cells and memory cells. The results that fol- 
low are averaged over 5000 instances of the ensemble. 

The average affinity of the population of antibodies 
evolves during mutation and selection. The evolution of 
the affinity energy under the two different strategies is 
shown in Fig. ^ The GSS+PM yields sequences with a 
lower energy than those from PM. The gap is due to each 
of the subdomains having an energy improved from an 
average to an exceptional one, AU « 10(C/^ st — ) = 
—0.8, and due to better sd-sd interactions found by GSS. 
In other words, the former is a better mechanism than 
the latter in searching sequence space for higher affin- 
ity antibodies to the given antigen. The best binding 
energy averaged over 5000 instances of the ensemble is 
U G = -21.9 in the GSS+PM case and U p = -19.7 in 
the PM only case at round 30, corresponding to affinities 
of K = 6.7 x 10 7 1/mol and K — 1.6 x 10 6 1/mol, re- 
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FIG. 2: Affinity of memory antibody sequences after a pri- 
mary immune response for the two different immune system 
strategies (PM and GSS+PM) to altered antigens. The bind- 
ing constant is K, and the antigenic distance of the new al- 
tered antigen from the original antigen is p. Cross-reactivity 
ceases at larger distances in the GSS+PM case (no cross- 
reactivity for p > 0.472) than in the PM only case (no cross- 
reactivity for p > 0.368). 



spectively. So, GSS+PM yields an affinity more than an 
order of magnitude improved, which is even better than 
the affinity obtained from a secondary response by PM 
only yj. That is, the multi-spin flips of GSS+PM espe- 
cially accelerate the optimization of u sd ~ sd in amino acid 
sequence space. Antibodies with a higher affinity to the 
antigen work more effectively in many ways. For exam- 
ple they can neutralize bacterial toxins, inhibit the infec- 
tiousness of viruses, and block the adhesion of bacteria to 
host cells at lower serum concentrations of antibody Q. 
Based solely on Fig. [2 it is hard to understand why Dar- 
winian evolution did not result in GSS+PM or any other 
more efficient strategy, rather than somatic hypermuta- 
tion, as the preferred strategy for B cell expansion. 

A calculation of cross-reactivity, which quantifies the 
specificity of the antibody, is used to compare the anti- 
bodies generated by the two strategies. Thirty rounds of 
primary response affinity maturation are conducted for 
both PM and GSS+PM. The antigen is then changed 
by p, and so each interaction parameter in the Hamilto- 
nian is changed with probability p. The affinity con- 
stants are calculated for the new antigen in the two 
cases. As shown in Fig. we find that cross-reactivity 
ceases approximately Ap — 0.10 later in the GSS+PM 
case. We also find that the affinity constant, K, de- 
creases exponentially with the degree of antigenic change, 
i.e. the binding energy, U, increases linearly with the 
change in the antigen, p. Within the region of specific 
binding (K > 100 1/mol), affinity is always better in 
the GSS+PM case. These cross-reactivity experiments 
show that the antibodies generated by the dynamics of 



GSS+PM can recognize more antigens and with higher 
affinity than those given by the PM dynamics only. Such 
cross reactivity has recently been observed 0. Inter- 
estingly, U(p = 1) = U(T = 0)/5, since at p = I the 
only correlation is that oti remains unchanged 20% of the 
time. The fraction p of new a values in u sd ~ sd and U c 
are centered around zero and negligible compared to the 
original, negative evolved values. In U sd , a fraction p of 
the Oi are changed, and the original, evolved sequence 
gives a value of roughly zero in a changed type ctj . The 
overall energy lost in U is, thus, proportional to p, Fig. [21 

Now the question we must answer is how many pro- 
tein molecules are there between pi = 0.368 and pi — 
0.472 that can be recognized by the antibodies produced 
through VDJ recombination and GSS+PM, but cannot 
be recognized by the antibodies produced through VDJ 
recombination and PM only. Antibodies recognize and 
bind to epitope regions of an antigen, so we are inter- 
ested in how many epitopes these two different classes 
of antibodies are likely to recognize. We now show that 
the abnormal antibodies recognize 10 3 times more epi- 
topes than do the antibodies produced by the normal 
primary immune response. We take the typical epitope 
length to be B = 20 amino acids 0, Q|. In the the- 
ory, the antigenic distance is defined by the chemical 
composition of the epitope. The total number of pos- 
sible epitopes is, therefore, 20 s since there are 20 differ- 
ent amino acids . In Fig. [31 we show the normal- 
ized number of epitopes that are at an antigenic dis- 
tance of p from the native epitope of the antibody, given 
by N(p) = WB\/[i\(B - where i = pB. This 
number is N( Pl )/20 B = 2 x 10~ 12 and N(p 2 ) /20 B = 
2 x 10~ 9 . The number of epitopes between p = 
and pi and the number between p\ and p 2 is approx- 
imately A\ w 2 x 10~ 12 and A 2 ~ 2 x 10~ 9 , respec- 
tively, since N grows exponentially with p. In the limit 
of large B, N{p 2 )/N{ Pl ) ~ exp{-B[ln(l-p 2 )/(l-p 1 ) + 
Pilnl9(l/pi-l)-p 2 lnl9(l/p 2 -l)]}. This ratio is 10 3 
for 20 amino acids. The ratio varies between 160 and 
4900 at the 15-25 amino acid limits of known epitope 
sizes. The ratio varies by only 30% with ±10% variation 
of both p values, the precision to which thegeneralized 
NK model p values agree with experiment |d \\(\ . 

In order to answer the question of why selection has 
favored the relatively slower dynamics in the immune re- 
sponse, we consider the function of the immune system 
as a whole rather than just the function of a single an- 
tibody binding to a single foreign antigen. The immune 
system is a delicate balance of a strong response to in- 
vading non-self antigens and weak or no response to self 
antigens. Thus, it is very important for the antibodies 
produced by the immune system to discriminate self from 
non-self |l7lllSlllflj . An immune system incapable of rec- 
ognizing an invading pathogen and initiating a response 
is an inadequate defense mechanism. On the other hand, 
production of antibodies binding to self antigens results 
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FIG. 3: The number of possible epitopes that are an antigenic 
distance of p from another epitope. The epitopes are assumed 
to be 20 amino acids in length. The plot is normalized, so that 
£ p iV( P ) = 20 20 . 

in autoimmune diseases, e.g. diseases such as type I di- 
abetes and rheumatoid arthritis. It might be thought 
that GSS+PM would be superior dynamics, if only the 
immune system were able to perform rs 10 rounds. How- 
ever, each subsequent exposure to related antigen would 
lead to another w 10 rounds. After a couple such ex- 
posures, GSS+PM has evolved an an unacceptably large 
binding constant, which saturates after several exposures 
to of the order of 10 10 1/mol, whereas PM saturates at the 
value of 10 7 1/mol. Even with PM dynamics, our model 
predicts that chronic infection may lead to autoimmune 
disease, a mechanism postulated to be responsible for 
some fraction of rheumatic diseases, including arthritis 
|20|. To resolve the significance of this mechanism, it 
would be interesting to see whether the distribution of 
onset times is broad, as predicted by our model. 

It has been experimentally shown that cross-reactivity 
can lead to autoimmune disease l\. From Figs. EH3 we 
know that the antibodies obtained at the end of a primary 
immune response can recognize a random epitope with a 
length of 20 amino acids with a probability of 10 -12 . 
It is also known that in a typical cell there are s» 10 4 
proteins, each with a length of as 500 amino acids |2lj . 
Antibodies binding to proteins recognize only surface epi- 
topes. The amino acids exposed on the surface of proteins 
are typically part of a loop or turn, and typically 1/3 of 
the loop or turn is exposed. Thus, a typical contiguous 
recognition fragment length is 6-7 amino acids. Given 
a typical length of 20 amino acids for the entire loop or 
turn, and the typical epitope size of 20 amino acids, an 
antibody will recognize roughly 3 non-contiguous regions 
of length 6-7 amino acids in the protein sequence of 500 
amino acids. There are approximately 500 3 ~ 10 8 such 
epitopes. Given the 10 4 proteins in a cell, there will be 
« 10 12 total epitopes expressed in each cell. Thus, the 
number of epitopes recognized by a typical antibody in 



each cell is A\ x 10 12 w 1. Even taking an exceptional 
protein of length 1000 amino acids, the number increases 
only to 10. 

This remarkable result shows that the antibodies pro- 
duced by the immune system recognize on average only 
their intended target. Conversely, antibodies that would 
be evolved through an immune response composed of 
VDJ recombination followed by a period of GSS+PM 
would recognize on average A^ x 10 12 » 10 3 epitopes in 
each cell. Such antibodies, while having higher affini- 
ties for the intended target, would also lead to many 
more instances of autoimmune disease. Such promiscu- 
ous antibodies would place too large a burden on the reg- 
ulatory mechanisms that eliminate the occasional aber- 
rant antibody j^. Thus, we find that selection has suc- 
cessfully evolved the human immune system to generate 
antibodies that recognize on average only the intended 
epitope after a humoral immune response. Inclusion of 
"more efficient" moves is generally excluded by the bound 
A x 10 12 = 0(1). 



[1] C. A. Janeway, P. Travers, M. Walport, and M. J. Shlom- 
chik, Immunobiology (Garland Publishing, New York, 
2004), 6th ed. 

[2] R. A. Goldsby, T. J. Kindt, B. A. Osborne, and J. Kuby, 

Immunobiology (W. H. Freeman and Company, New 

York, 2002), 5th ed. 
[3] K. H. Fischer and J. A. Hertz, Spin Glasses (Cambridge 

University Press, New York, 1991). 
[4] D. J. Gross and M. Mezard, Nucl. Phys. B 240, 431 

(1984). 

[5] M. Ney-Nifle and A. P. Young, J. Phys. A 30, 5311 
(1997). 

[6] V. Azcoiti, E. Follana, and F. Ritort, J. Phys. A 28, 3863 
(1995). 

[7] D. J. Earl and M. W. Deem, Proc. Natl. Acad. Sci. USA 

101, 11531 (2004). 
[8] L. D. Bogarad and M. W. Deem, Proc. Natl. Acad. Sci. 

USA 96, 2591 (1999). 
[9] M. W. Deem and H.-Y. Lee, Phys. Rev. Lett. 91, 068101 

(2003). 

[10] S. Kauffman and S. Levin, J. Theor. Biol. 128, 11 (1987). 
[11] R. S. et al, J. Mol. Bio. 263, 551 (1996). 
[12] M. G. Kidwell and D. R. Lisch, Evolution 55, 1 (2001). 
[13] P. D. Holler, L. K. Chlewicki, and D. M. Kranz, Nature 

Im. 4, 55 (2003). 
[14] E. T. Muiioz and M. W. Deem, Vaccine 23, 1142 (2004). 
[15] N. A. Campbell and J. B. Reece, Biology (Benjamin 

Cummings, San Francisco, 2004), 7th ed. 
[16] V. Gupta, D. J. Earl, and M. W. Deem (2005), q- 

bio.BM/0503030. 
[17] J. K. Percus, O. E. Percus, and A. S. Perelson, Proc. 

Natl. Acad. Sci. USA 90, 1691 (1993). 
[18] V. Detours, H. Bersini, J. Stewart, and F. Varela, J. 

Theor. Biol. 170, 401 (1994). 
[19] J. Sulzer, L. van Hemmen, A. U. Newmann, and U. Behn, 

Bull. Math. Biol. 55, 1133 (1993). 
[20] M. Leirisalo-Repo, Curr. Opin. in Rheumatology 17, 433 



■5 



(2005). (2005). 
[21] T. Tan and M. W. Deem, Physica A 350, 52 (2004). 
[22] M. Larche and D. C. Wraith, Nature Med. 11, S69 



